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Abstract 

We investigate the time evolution of the heteropolymer model in- 
troduced by Iori, Marinari and Parisi to describe some of the features 
■ of protein folding mechanisms. We study how the (folded) shape of the 

chain evolves in time. We find that for short times the mean square 
distance (squared) between chain configurations evolves according to 
a power law, D ~ t v . We discuss the influence of the quenched dis- 
order (represented by the randomness of the coupling constants in 
the Lennard- Jones potential) on value of the critical exponent. We 
find that v decreases from | to \ when the strength of the quenched 
disorder increases. 
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1 Introduction 



Modeling the protein folding is one of the most relevant challenges open to 
Statistical Mechanics. The idea that relevant issues of the folding process can 
rely on the disordered nature of the amino-acid system has been put forward 
and investigated by many groups (see for example refs. [0J and references 
therein). Shakhnovic and Gutin@ in their seminal work have been giving a 
first quantitative form to these ideas, implementing the mechanism of Parisi's 
Replica Symmetry BreakingS. Mezard and Parisi have applied a systematic 
RSB treatment to random manifolds^. 

Iori, Marinari and Parisi in ref. || (IMP) have studied a very simple 
heteropolymer model, attempting to emphasize the difference between the 
static structure of spin glasses and native proteins. The model turned out 
to present potential important features, like the dominance of one or a small 
number of ground states. In ref. || the model is studied in d = 2, and ref. 
[0] studies the case of Lennard- Jones homopolymers. 

The Hamiltonian of the IMP model (describing a self-interacting het- 
eropolymer) has the form 

where i,j range from 1 to N (the number of sites of the polymeric chain), 
and d?j = (r$ — r,) 2 . The first term represents a harmonic force holding 
the adjacent sites together. This term prevents the chain from breaking into 
pieces. The contributions proportional to A and R (attractive and repul- 
sive) form a conventional Lennard- Jones (L-J) potential. The last term is 
a quenched disorder contribution; it models a nonuniformity of the dipoles 
interacting by the L-J potential. The values of rjij are independent and un- 
corrected stochastic variables (which do not vary in time). The rj^/s have 
mean value zero and variance one. The set of coefficients rjij can be viewed 
as characterizing a single native protein. 

The most relevant questions concern typical realizations of the heteropoly- 
meric chains. In the following we will mainly interested in quantities which 
are averaged both on the thermal noise and on different realizations of the 
T)ij couplings. 

In this note we investigate the short time behavior of the heteropolymeric 
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chains whose dynamics is defined by the Hamiltonian (|I|). Specifically, we 
study the diffusion of the distance between two chain configurations defined 

as 

a^/^^E^-^?) 2 , (2) 

where a and (3 label two chain configurations at different instances of 
time. The time evolution of D 2 tells us about the diffusion of the shape (and 
size) of the folded chain. We concentrate our attention on the dependence of 
diffusion on e, the quenched noise strength. 

Analytical studies of the diffusion in the (simple) case of polymer where 
the interaction is local along the chain are possible in the context of normal 
mode decomposition of the Langevin dynamics (LD). Such simple models 
are known to yield a characteristic _power law dependence of the mean-square 
displacement (squared) over timeO, £2. This power law holds in the range 
of times where the finite size effects are still invisible, but collective effects 
are already important®. This is an interesting regime to which we will 
dedicate attention in the following. Other relevant models are discussed 
in the context of the problem of surface aggregation (see for example || 
and |I0[). For example for a nonlinear model with a local interaction in 

2 

'1 + 1) dimensions the mean square displacement (squared) grows as is in 



the relevant regime!^ Hi) 

In the case of interest (heteropolymers with long range interaction on the 
chain, and a L-J potential with a random quenched attractive contribution) 
the problem is very complex (since the Hamiltonian is disordered and non- 
local, and frustration can play a role). It is also important to remind that we 
want to study a situation where N is large but not infinite, and this finiteness 
can play a role. This is true both for our model (we will present here results 
obtained with iV = 30) and for the potential application to protein folding. 
In a biological application of our model the sites of our chain have to be iden- 
tified with (pre-assembled) parts of the secondary structure, and N would be 
of order 10. The heteropolymeric chain is finite and localized (at equilibrium 
the distance between chain sites is finite: translational and rotational degrees 
of freedom do not enter in our definition of mean displacement). This im- 
plies that, after a time large enough, the mean square displacement reaches a 
maximum (which, roughly speaking, characterizes the amplitude of thermal 
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fluctuations and depends on spatial extent of the chain). 

Nonlinear systems with external bias are studied numerically and are 
expected to yield crossover between t 1//2 and t 2 ^ 3 power laws. It is not com- 
pletely clear whether the apparent nonlocality of interaction can be treated 
effectively as a bias (an effective mean field due to cloud of particles) incor- 
porated in a local nonlinear model describing nearest-neighbor interaction 
along the chain. 

2 The Equilibrium Shape 

We have tried to get some more informations about the equilibrium shape 
of the IMP heteropolymers. We have mainly used in our runs the same 
parameters used in imp!. 

Let us look at some raw number in order to try to understand what 
happens when we go from the homopolymer phase to the folded phases, 
where the disorder plays a crucial role. In absence of the quenched disorder 
(e = 0), for {3 = 1, N = 30, a = 3.8, r = 2 and h = 1, the total chain 
energy is < Et > = —273.0 ± .5, the average chain first neighbor square 
length < d 2 i+1 > = 1.55 ± .10 , and the end to end chain length < d 2 N1 > 
= 6.0 ± 1.5. The corresponding values for e = 6 are < Ej- >= —635 ± 50, 
< d 2 i+1 > — 1.7 ± 0.5 <d 2 Nl > = 3.5 ± 3.5 (here the errors are computed by 
averaging over different realizations of the rjij couplings). These number are 
here only to hint an order of magnitude. Indeed in ref. || it has been seen 
that 100 millions of full chain sweeps are not at all sufficient to explore the 
entire phase space. 

The values of e needed in IMP (for iV = 30) to go deep in the folded 
phase (e ~ 6) are quite large. We have then to be quite careful in checking 
which terms are contributing to the ground state total energy. We find that 
in both the disordered and in the e = regimes the attractive part of the 
L-J energy expectation value (in which we include the part proportional to 
e) is roughly equal to twice the repulsive part of the energy. This is what we 
find for the minimum of a L-J potential for a two-particle system, where at 
equilibrium 



(d is the distance between the two particles at equilibrium). For these 
values of the parameters the harmonic term contributes about ten percent 
to < E T > but practically does not affects the value of d. 

For e = in the limit N — > oo the ground state configuration is such 
that the chain sites lie on a face-centered cubic lattice. The lattice (which 
in this case coincides with the solution to the kissing problem, i.e. finding 
the maximum number of spheres of diameter d tangent to a given one) can 
be seen as made of centers of spheres of diameter d packed regularly with 
the maximum density (which in 3d is 12). In this limit the harmonic energy 
term only determines the order in which the chain sites are placed on the 
lattice sites. Fukugita, Lancaster and Mitchard0 have shown that in d = 2 
the lattice structure is quite clear in the e = case, and survives (although 
the evidence of ref. || is in this case less compelling) the transition to the 
disordered phase. In the 3c? case we find that for N = 30, in the folded 
phase, the lattice structure is lost. We clearly see the lattice structure in 
the e = case, and we see that it becomes weaker but partly survives when 
thermalizing the chain in the strongly disordered potential ||12|| . 



3 The Time Dependence of the Shape 

In order to discuss the diffusion of chain configuration quantitatively, it is 
more convenient to introduce the following definition of distance between two 
configurations 

^(a^^E^-^) 2 , (4) 

Since 

(4f -d^) = (d%-d^)(d% + d^), (5) 

where the second factor is much larger than the first, the exponent char- 
acterizing the power law short-time dependence of D2 and -D4 is the same. 
One could also use the two definitions of the distance introduced in ref. [|], 
and expect to find the same critical behavior. We study the diffusion of the 
polymer shape by analyzing an ensemble average of -D4 as a function of the 
time separation between the two configurations. We average over the time 
dynamics and over different realizations of the 77^ couplings, and compute 
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D 4 (t) =<D 4 (a(t),a(0))> , (6) 

where a(0) is an equilibrium configuration. 

In all our runs we have used a standard Monte Carlo dynamics. 

The fluctuations of the critical exponent governing the behavior of D 4 (t) 
from sample to sample are quite small. We average D 4 over a few hundreds 
(100-500) of starting points for a given sample. 

In order to check that we are really picking up the correct universal be- 
havior we have studied a simple harmonic chain. Let us stress that the 
problem is very delicate: there is a non universal short time region (where 
the discreteness of the Metropolis procedure plays an important role), an 
asymptotic constant value for D 2 (with some universal approach to D 4 , and 
in between the (short) time region we are interested in. We have to check 
that we are observing, in this region, a true scaling behavior, and we can be 
sure of that only in the limit of large N. In order to check our results for the 
simple harmonic chain we have compared them to the analytic expression 
resulting from the normal mode expansion of the Langevin dynamics (LD). 

The Langevin equation for the harmonic chain is of the form 

(IT ■ 

= -k {2n - r l+1 - ri _i) + ft , (7) 
where the uncorrelated random forces fj satisfy 

< fMfjih) >=2(k B T 6(t! - t 2 ) 5 itj . (8) 

(, a friction constant, sets here the time scale for the problem. In the 
continuum limit (2r^ — r i+1 — rj_i) may be replaced by d 2 r/dn 2 . One can 
introduce normal modes in the standard way by 




Using Wick's theorem for evaluating average of products of the coordi- 
nates <XjXkX m x n > one can easily compute the D^t), finding 

D4(t) = £K,(o) - iUj(t)] , (10) 

where 
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= c {£ e ^l(cos(^) - cos(^)) 2 } 2 , (11) 

p // 1 V iV 

and 



64x4^ = 64x^1; (12) 
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iV 6 C^ 2 



" A;T3vr 2 /ivr 2 ' (13) 

b 2 = kT/3h is the asymptotic equilibrium expectation value of d 2 i+l . We 
show the comparison of our numerical fitted data with the analytic result in 
fig. 1. We can see very well that the power for short times is \. Saturation 
starts to show up for larger times (and it is very clear for very large times, 
not included in the figure), also if the statistical error is already becoming, 
in this regime, very large. 

Let us just remind again that we are using a discontinuous dynamics, the 
Metropolis algorithm. We know that we are getting the correct asymptotic 
equilibrium distribution, but for short times we do not have a simple corre- 
spondence with the corresponding Langevin dynamics. Let us discuss this 
point in some detail. Consider a Metropolis dynamics, where we propose 

— » 

a trial random increment defined by a displacement vector 5 chosen from 
some probability distribution P(5). The increment is accepted with prob- 
ability P a = mm(l,eP( H ^~ H ( r+s ^). Therefore the average of a single step 
displacement is 

<r'-r> ~ <5 min(l,e /3{H{?) - H(?+ ^ )) )> . (14) 

For small (35 ■ Vif(r), (corresponding to a large acceptance factor) the 
< r i — f> b e approximated as 

-(3 <5B(5- VH(r)) (5- VH(r))> = ~ <5 2 >VH(r) , (15) 

where 6 is the step function. This is the same form one gets for a Langevin 
step, with a scale factor which depends on the acceptance ratio of the Metropo- 
lis procedure, < 5 2 > (3. Normally one uses the Metropolis algorithm with 
a trial displacement that is not small: this is indeed the big advantage of 
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the Monte Carlo method. That means that in general condition ( |T5| ) does 
not hold. We expect the Metropolis dynamics to reproduce the continuous 
dynamics only for times larger than the typical time of continuous dynamics 
necessary to reach the mean displacement < 5 >. This is the time region 
we are interested in, and that we discriminate before analyzing the dynam- 
ical critical exponent. In this region (of times large enough not to feel the 
non-universal details of the dynamics, but small enough not to describe the 
relaxation to the asymptotic mean displacement) we will see that we can 
determine a critical exponent v 

D 4 (t)~t u . (16) 

In the context of Langevin Dynamics we can relate the exponent v with 
the asymptotic behavior of dynamics linearized around local minima. We 
consider the normal modes eigenvalues of the linearized dynamics, X p . Let 
us assume that asymptotically (i.e. in the continuum limit of large N and 
small p) 

\ = V a (17) 

and that the fluctuations of all modes are the same and not correlated. 

The width of the probability distribution of the p-th mode behaves as — §=, 

V -V 

where C is a mode independent constant. Then we obtain 

D*{t) = EKi(°) - > ( 18 ) 

where 

p V 

where the the scalar products of the p-th mode vector x p with 

the position vector of the i-th site, r^. Changing the order of the summation 
and converting the sum over p into an integral (in the continuum limit) we 
obtain 

D 4 (t) ~ / dp p- a (l - e~ p ) ~ t 1 -- , (20) 
Jo 
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where the coefficients c P) ; were replaced by their average values and summed 
over |J. Therefore, under the above approximation, the expected relation 
between the exponent a% and the exponent v is 



1 



(21) 



v = 1 



a: 



We have verified numerically that in the large e region the above relation 
is actually satisfied. This is consistent with the fact that v is a true universal 
exponent, dependent only on the energy spectrum and independent from j3. 

Our numerical result do clearly exhibit, in the large e region, the scaling 
relation 



for N as small as 15— 30, justifying the claim we are close to the continuum 
limit (obviously we ignore the first six eigenvalues which are zero due to the 
rotational and translational symmetries). 

4 Numerical Estimates of v 

We present here our numerical determination of the dynamical critical ex- 
ponent v we have defined before, from the short (but not so short) time 
behavior of D^. 

The four parameter which determine the Hamiltonian, A, R, h and e 
affect the behavior of D^(t). Let us note, at first, that the time scale r r after 
which the power law does not hold any more 



is a decreasing function of the strength of interaction, and an increasing 
function of the system size. Secondly, the coefficient k decreases with the 
strength of coupling and increases with the temperature and the system size. 
In fig. 1 we have shown the behavior of an harmonic chain, together with 
the results of the theoretical analysis. 

In fig. 2 we show the time behavior of ordered systems for two different 
values of the parameters. Here (in the globular phase) we find a good fit with 
v — |, which coincides with the prediction of models with local interaction. 




(22) 



D 4 (t) ~kt u ,t « T r 



(23) 
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The fact that this exponent stays the same even if we increase the local 
interaction be increasing h ten times while keeping R and A fixed suggests 
that in the coil phase the diffusion properties in the absence of the quenched 
noise can be adequately described by a model of sites with a local interaction 
placed in an effective external field. 

In the following we will mainly be interested in the dependence of the 
exponent v over e. 

Let us note at first that the exponent v does not seem to depend on the 
particular realization of the quenched noise. We have not performed large 
scale simulations, which would be needed in order to get quantitative precise 
results in systems with such a dramatic critical slowing down (see ref. ||), 
but we have got from our simulations quite a precise qualitative picture. We 
show in fig. 3 D±(t) for three different realizations of the noise, with a very 
similar power behavior. The non-universal coefficient of the power behavior 
does indeed depend on the given noise realization, but v does not. Of course, 
the energy of interaction itself does not determine v (as it does not determine 
the equilibrium shape). 

In fig. 4 we show D±(t) for different e values. For increasing e we observe 
a reduction of v: for e = 10 v — .52 ± .05 

As noted in ref. || the correlation time shows a dramatic dependence 
on e: in fig. 5 we show the distribution probability Pt{D$) for four values 
of e. We have selected the same time moment for the 4 distributions. The 
probability PtiD^) is obtained by averaging over several hundred initial con- 
figurations. While for e = the P{D/±) for times ~ 100 ( where 9 = 10 4 
updates of each chain site) is nearly stationary, for e = 6 times of two orders 
of magnitude longer are may be just starting to be enough (see fig. 6). The 
new Tempering approach to Monte Carlo dynamics, recently proposed in ref. 
[[Uj, has chances to alleviate the situation. 
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Figure Captions 



1 Comparison between the D±(t) (here tuc is the number of full MC sweeps 

of the chain) estimated in our Monte Carlo runs for the purely harmonic 
chain (scattered dots) and the theoretical prediction from the Langevin 
dynamics normal mode expansion (Eq. pi], continuous curve). Here 
T r = 800. 

2 As in fig. 1, but for a homopolymer (e = 0) with Lennard- Jones inter- 

action, with A = 3.8 and R = 2, log-log plot. The two curves are for 
h = 1 and for h = 10. 

3 -D 4 (t) for three different realizations of the quenched noise rjij, e = 6. 

a) linear plot, b) log-log plot (arbitrary normalization). 

4 \og(D 2 ) versus log(t) for different values of e (arbitrary normalization). 

The steepest line (labelled with random) corresponds to a system of 
free particles - all proposed moves accepted. 

5 The distribution probability Pt{D^) for four values of e. We take the same 

time for the 4 distributions (twice the maximum time shown in fig. 4). 

6 The plot of time dependent probability distributions for e = 6 for three 

different values of time. The smallest time t corresponds to the time of 
fig. 5. 
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